

Compressed Semi-Discrete Central-Upwind Schemes for 
Hamilton-Jacobi Equations 

Steve Bryson*, Alexander Kurganov*, Doron Levy* and Guergana Petrova^ 


Abstract 

We introduce a new family of Godunov-type semi-discrete central schemes for multidimensional 
Hamilton-Jacobi equations. These schemes are a less dissipative generalization of the central- upwind 
schemes that have been recently proposed in series of works. We provide the details of the new family 
of methods in one, two, and three space dimensions, and then verify their expected low-dissipative 
property in a variety of examples. 

AMS subject classification: Primary 65M06; Secondary 35L65 

Key Words: Hamilton-Jacobi equations, central-upwind schemes, semi-discrete methods. 


1 Introduction 

We consider the multidimensional Hamilton-Jacobi equation, 

<Pt + H{V x tp) = 0 , x € (1.1) 

with Hamiltonian H . First-order numerical schemes that converge to the viscosity solution of (1.1) 
were first introduced by Crandall and Lions in [5] and by Souganidis in [17]. Recent attempts to 
obtain higher-order approximate solutions of (1.1) include upwind methods, discontinuous Galerkin 
methods, and others. Here, we study a class of projection-evolution methods, called Godunov-type 
schemes. The main structure of these schemes is as follows: one starts with the point values of the 
solution, constructs an (essentially) non-oscillatory continuous piecewise polynomial interpolant, and 
then evolves it to the next time level while projecting the solution back onto the computational grid. 
The key idea in Godunov-type central schemes is to avoid solving (generalized) Riemann problems, by 
evolving (locally) smooth parts of the solution. 

Second-order staggered Godunov- type central schemes were introduced by Lin and Tadmor in [14, 15]. 
L l -convergence results for these schemes were obtained in [14]. More efficient non-st aggered central 
schemes as well as genuinely multidimensional generalizations of the schemes in [15] were presented in 
[1], with high-order extensions (up to fifth-order) proposed in [2, 3]. 
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Second-order semi-discrete Godunov-type central schemes were introduced in [12], where local speeds 
of propagation were employed to reduce the numerical dissipation. The numerical viscosity was fur- 
ther reduced in the central-upwind schemes [10] by utilizing one-sided estimates of the local speeds 
of propagation. Higher-order extensions of these schemes were introduced in [4], where weighted es- 
sentially non-oscillatory (WENO) interpolants were used to increase accuracy. WENO interpolants 
were originally developed for numerical methods for hyperbolic conservation laws [16, 8], and were first 
implemented in the context of upwind schemes for Hamilton- Jacobi equations in [7]. 

Godunov-type central-upwind schemes are constructed in two steps. First, the solution is evolved to 
the next time level on a nonuniform grid (the location of the grid points depends on the local speeds, 
and thus can vary at every time step). The solution is then projected back onto the original grid. The 
projection step requires an additional piecewise polynomial reconstruction over the nonuniform grid. 
In this paper we show that in the semi-discrete setting different choices of such a reconstruction lead 
to different numerical Hamiltonians, and thus to different schemes. In particular, we can recover the 
scheme from [10]. A more careful selection of the reconstruction results in a new central-upwind scheme 
with smaller numerical dissipation. This approach was originally proposed in [11], where it was applied 
to one-dimensional (1-D) systems of hyperbolic conservation laws. It has been recently generalized and 
implemented for multidimensional systems of hyperbolic conservation laws in [9]. 

The paper is organized as follows. In §2, we develop new semi-discrete central-upwind schemes for 
1-D Hamilton- Jacobi equations. We also review the interpolants that are required to complete the 
construction of the second- and fifth-order schemes. Generalizations to more than one space dimensions 
(with special emphasis on the two-dimensional setup) are then presented in §3, where the corresponding 
multidimensional interpolants are also discussed. In §4, we evaluate the performance of the new schemes 
with a series of numerical tests. Finally, in the Appendix, we prove the monotonicity of the new 
numerical Hamiltonian. 

2 One-Dimensional Schemes 


2.1 Semi-Discrete Central-Upwind Schemes for Hamilton- Jacobi Equations 

In this section, we describe the derivation of a new family of semi-discrete central-upwind schemes for 
the 1-D Hamilton- Jacobi equation, 

<pt 4* H i}px) = 0, x € R, (2*1) 

subject to the initial data ip(x,t = 0) = ipo{x). We follow the approach in [10] (see also [12]). For 
simplicity we assume a uniform grid in space and time with grid spacing Ax and At, respectively. The 
grid points are denoted by Xj := j Ax, t n := nAt, and the approximate value of ip ( Xj } t n ) is denoted by 

Assume that the approximate solution at time t n , <p”, is given, and that a continuous piecewise- 
polynomial interpolant <p(x, t n ) is reconstructed from tp™. At every grid point, the maximal right and 
left speeds of propagation, a t and aj , are then estimated by 


a+ = _ max _ (iT(u), 0} , a j 

minlv?! wt}<u<max.{ip x ,<£>+} 


min { H \u), 0 } , 

min {(p x ><pt}<u<max.{<p x ,<pt} 


( 2 . 2 ) 


where tp^ are the one-sided derivatives at x = Xj % that is 


<Px '■= <Px{Xj±Q,t n )- 
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Figure 2.1: Central-upwind differencing: 1-D 


If the Hamiltonian is convex, (2.2) reduces to 

a t ~ max {H'(<Px)>H'( l Px)>Q} > a" = |min{ir , (v7*),/Z' / ( ¥ ?+) ) 0}|. (2.3) 

We then proceed by evolving the reconstruction Ip at the evolution points x™ ± := Xj ±afAt, to the next 
time level according to (2.1). The time step At is chosen so that for all ;. Therefore the 

solution remains smooth at x^ ± for t G [ t n , t n+1 ] (see Figure 2.1) and we can compute the values of the 
evolved solution {p^ 1 } by the Taylor expansion: 

v # 1 = £(*?±, n - A tH (v x (x] ± , n) + O (At 2 ) . (2.4) 

Using the values on the nonuniform grid {a;™ ± }, we construct a new quadratic interpolant 

^(£,£ n+1 ) on the interval 






x’-j_ — X'/ 
3 + 3 - 


(x - x^_) 4- -((p x X )] +I (x - x”_)(x - x^ + ), 


where (£ IS )” +1 is yet to be determined and is an approximation to p X x{%ji £ n+1 ), = (x™ + + x^_)/2. 

The projection back onto the original grid is then carried out by evaluating ^(x,^ 4 " 1 ) at Xj, 

<P? +1 ■■= *" +1 ) = - l(?xi )” +1 ^aT (At) 2 . (2.6) 


a j +a j 


a j + a j 


Note that if the Riemann fan is symmetric, that is, if at = , then X™ = Xj. Substituting (2.4) in 

(2.6) yields 


°? +1 = -^r(ip(x^,t n )-A tH(v x (x?_,t n ))) 

+ -r^-=(<p(x] + ,n - AtH(ip x (x] + ,n)) - l(^), n+ 1 at a -(At ) 2 + <D(At? 


Using the Taylor expansion, 

W±, t n ) = ± Ataftfi + O(At ) 2 


( 2 . 8 ) 
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we arrive at 


<P 


TL -{-1 


= <P? + A ('fit - <PZ) ~ -+^-r \ a J H(<p x (Xj+> t n )) + afH(lp x (x?_,t n ))] 

Qj 2 “r CL- CL- CL- L J 


a i + ~i 


--(tp xx )] +l ataJ(Aty + 0(Aty. 

We then let At — ► 0, and end up with the (family of) semi-discrete central-upwind schemes: 


d a j H(<p+) + afH{<p x ) _ / - tp~ 1 

— ipi(t) = — - z x + aja , Af — xr - x lim 

dt dj ■+- dj \ &j + dj 2 At— *0 




(2.9) 


( 2 . 10 ) 


Here, the one-sided speeds of propagation, a are given by (2.2), and are the left and right derivatives 
at the point x = Xj of the reconstruction £>(*,£) at time t. 

Finally, in order to complete the construction of the scheme, we must determine (^p X x)] +l 
example, selecting {tpxxYj^ 1 t 0 independent of At gives 


For 


A lim Q [A t(<p xx )y +1 = 0, 


and then (2.10) recovers the central-upwind scheme in [10]. However, since the interpolant ^(*,£ n+1 ) 
is defined on the intervals whose size is proportional to At, it is natural to choose ({ p xx )j +1 

to be proportional to 1/At. In this case, the approximation of the second derivative in (2.10) will 
add a non-zero contribution to the limit as At — > 0. At the same time, to guarantee a non-oscillatory 
reconstruction, we should use a nonlinear limiter. For example, one can use the minmod limiter: 


At((p xx )^ 1 = 2minmod 


'fix (z" + ,i n+1 ) - V>i(s?,f t+1 ) V> x (x?,t n+1 ) - PxjXj-, t n+1 ) ' 


CL^~ -f- CL ■ 
3 3 


a+ + aj 


( 2 . 11 ) 


where is the derivative of (2.5), tp x (x™ ±% t n+1 ) are the values of the derivative of the evolved recon- 
struction <£>(•, £ n ) at t = t n+1 , and the multivariate minmod function is defined by 


min{xj}, if Xj > 0 V?, 

3 

minmod(xi,X 2 , ...) := { max{xj}, if Xj < 0 Vj, 

i j 

0, otherwise. 


( 2 . 12 ) 


A different choice of limiter will result in a different scheme from the same family of central-upwind 
schemes. 

All that remains is to determine the quantities used in (2.11). Since all data are smooth along the 
line segments (x^ i5 £), t n < t < £ n+1 , we can use a Taylor expansion to obtain 

ip x (x^ ± ,t n+1 ) - ip x (Xj ± ,t n ) + O(At). 

According to (2.5), the derivative ip x (x^, £ n+1 ) of the new reconstruction ip at time level f n+1 is 


(2.13) 


P n t l ~ ¥>" +1 
J ( aj + dj )At 


and after substituting (2.4) and (2.8) into (2.14), we obtain 

'Tx *- L \YZ\*j+l t ' 

(dj + dj ) (at + dj ) 


Mxj,t n+1 ) - _ - H '(^( j "+. tn ))- g (^(^_^ n )) 


+ 0{At). 


(2.14) 


(2.15) 


i 
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Passing to the semi-discrete limit (At — > 0) in (2.11), (2.13), and (2.15) gives 


lim 

At->0 




n-f I 


(at + dj ) 


— minmod (ip+ - ^ nt , - <p x ) , 


where 


:= lim [Mx?,t n+1 )} = 

J (at + a j ) (aj + a j ) 


(2.16) 

(2.17) 


Finally, substituting (2.16) into (2.10), we obtain the 1-D low- dissipative semi-discrete central-upwind 
scheme: 


“ a j H (vi)+aj H (<p x ) , +_ 

<Pj(t) = — n + at a. 


d_ 

dt 


at + a j 


3 3 


<Pt 


’<Px 


at 4- a • 
3 3 


— minmod 


V at + aj 


,int 

X 


a t + aj 


> (2-18) 


where ip^ is given by (2.17). For future reference, we denote the RHS of (2.18) by —H BKLP . 

Notice, that in the fully-discrete setting the use of the intermediate quadratic reconstruction £ n+1 ) 
at level £ n+1 (as opposed to the intermediate piecewise linear reconstruction in [10]) increases the 
accuracy of the resulting fully-discrete scheme: <D({At) 2 + (Ax) r ) versus 0(At + (Ax) r ), where r is the 
(formal) order of accuracy of the continuous piecewise polynomial reconstruction £?(•, t n ). When we pass 
to the semi-discrete limit (At — > 0), both quadratic and linear interpolation errors go to 0, and therefore 
the (formal) order of accuracy of both (2.17)-(2.18) and the semi-discrete scheme in [10] is 0((A:r) r ), 
and the temporal error is determined solely by the (formal) order of accuracy of the ODE solver used 
to integrate (2.18). However, the minmod limiter introduces a new term that leads to a reduction of 
the numerical dissipation without affecting the accuracy of the scheme. To demonstrate this, we show 
that is always in the interval [min{^+, <£>“}, max{<p+, }], and therefore the absolute value of the 
term (<p+ — <^“) in the numerical dissipation in the scheme from [10] is always greater than the absolute 
value of the new term, that is 


\vi - v>x | > V>i - ¥>X - minmod (ip* - V4 nt > Vd nt - <P X ) | • 

Indeed, we have 

<int _ Vx + Oj <Px _ HQpj) - Hif£f _ ) _ . 

1 (af + aj) (af + aj) 


’af-H'iS) 



. a j+ a J . 

+ <Px 

. a t + a J . 


(2.19) 


where £ € (min{</>+, <p x }, ma x{tpt, Tx })• It follows from the definition of the local speeds (2.2) that 


at - J?'(0 > o, H\0 + a- > 0. 

Thus, (2.19) is a convex combination of ipf and y>“, and therefore ^ nt e [min{<^+,<pj},max{<^+,^-}] 
It was shown in [4] that the numerical Hamiltonian H KNP from [10] is monotone, provided that the 
Hamiltonian H is convex. Here, we state a theorem about the monotonicity of H BKLP — the new, less 
dissipative Hamiltonian in (2.18). The proof is left to the Appendix. We will consider only Hamiltonians 
for which H' changes sign, because otherwise either a~ = 0 or a + = 0 and the Hamiltonian in (2.18) 
reduces to the upwind one for which such a theorem is known. 

Theorem 2.1 Let the Hamiltonian H £ C 2 be convex and satisfy the following two assumptions: 

(Al) The function 

G(u, v) := 2 H"(u) [{u - v)H'{v) - (. H(u ) - £T(t>))] + [H'(u) - H'(v)\ 2 < 0 (2.20) 
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for all u and v in the set S ( u,v ) U S + (u, v), where 


, v) := | {u, v) : 


S ( u , 


:= < (ifc.v) : 


H(u) - H(y) H'(u ) + H\v) 


u — V 


}' 


< ~~ x ' -/ ; " x u > vf > v 


u — v 


( 2 . 21 ) 


and u* is the only point such that H f (u*) = 0; 

(A 2) For any v and for an arbitrary interval [a, b\, the sets S~(u, v) n [a } b] and S + (u,v) D [ a,b ] are 
either the empty set or finite unions of closed intervals and/or points . 

Then the numerical Hamiltonian in (2.18): 

( 2 . 22 ) 

where 

■ int a + u + + a~u~~ H(u + ) - H(u~) 

U := (a++o-) (a+ + a~) ’ 

and a + a + (u + ,u~) = max{H'(u+) } H f (u~~) y 0}, a” := a~(u + ,u”) = | min{if / (tt + ), 0}| is 

monotone, that is, H BKLP is a non-increasing function of u + and a non- decreasing function of u~. 


a + a 


— u ,( u 4 * — u int u int 

minrnod - 


u 


u 4 " -f- a 


y a^~ -f- a a~^~ -f* a 


Remarks. 

1. The classification of all Hamiltonians that satisfy conditions (2.20)-(2.21) is an open problem. 
However, we were able to verify these conditions for certain Hamiltonians. For example, a 
straightforward computation shows that for any convex quadratic Hamiltonians H(u) = au 2 + bu 4- c, 
the function G(u,v) = 0, the sets in (A 2) axe either 0 or one closed interval, or one point, and 
therefore the theorem holds. 

Another example, for which Theorem 2.1 is valid, is H(u) = u 4 . In this case, the sets (2.21) are 
S~(u y v) = {(u, u) : u + v > 0 > v} , S*(u,v) = {(u,v) : u + v < 0 < v} , 

and, as one can easily verify, the function 

G(u , v) = — 8(u — v) 3 (u 3 4- 3 u 2 v -f 6uv 2 -F 2v 3 ) < 0, 

in S~(u y v) U S p (u,v). As for the sets in assumption ( A2 ), they are either 0 or one closed interval, or 
one point. 

2. Notice that assumption (A 2) in Theorem 2.1 is needed only for technical purposes and in fact it is 
satisfied by (almost) every Hamiltonian H that arises in applications. 


2.2 A Second-Order Scheme 

A non-oscillatory second-order scheme can be obtained if one uses a non-oscillatory continuous piecewise 
quadratic interpolant ip. The values of the one-sided derivatives of ip at ( Xj,t n ) in (2.17) and (2.18), 
are given by 


= 


Ax 


Ax 


, \ m 

T "~T~ (.Tvxjj+i j 


( A P)"+I : = ¥>”+i “ ( P]> 


(2.23) 
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where the second derivative is computed with a nonlinear limiter. For example, 

(A^ + 3 - (A ?)» +1 (A^ + 3 - (A<p)"_i 


(<Ar 


'j+ 


— minmod 0- 


(Ax) 2 


2(Ax) 2 


(Av)" x-(A^_r 

J^2 J 2 

(Ax) 2 


.(2.24) 


Here, 0 E [1,2] and the minmod function is given by (2.12). The scheme requires a second-order ODE 
solver. 


2.3 Higher- Order Schemes 

In this section, we briefly describe the third- and fifth-order weighted essentially non-oscillatory (WENO) 
reconstructions. They were derived in [4] in the context of central-upwind schemes, and are similar to 
those used in high-order upwind schemes [7]. 

In smooth regions, the WENO reconstructions use a convex combination of multiple overlapping 
reconstructions to attain high order accuracy. In nonsmooth regions, a smoothness measure is em- 
ployed to increase the weight of the least oscillatory reconstruction. Here, we reconstruct the one-sided 
derivatives j at x — Xj for k = 1, . . . , d stencils, and write the convex combination: 


<Px 


Y w k,j(d 




X>t = 1 > 


w 


k,j 


>0, 


k = 1 


k=l 


(2.25) 


where the values are to be used in the scheme (2.17)~(2.18). The weights ■ are defined as 


a 


w\ 


k,j 


kj 


d 

I 

z=i 


E4 



(2.26) 


The constants are set so that the convex combination in (2.25) is of the maximal possible order of 
accuracy in smooth regions. We take p = 2 and choose e — 10 - * 6 to prevent the denominator in (2.26) 
from vanishing. 

A third-order WENO reconstruction is obtained in the case d = 2 with 


= 


<Pj+i - <Pj- 1 , _•> _ -<Pj - 2 + 4^j_i - 3 tfij 

2 Ax ’ 2 Ax 

-3 fj + 4 <pj+i - <£j+2 , _ <Pj + 1 - <£j-i 

2Ax ’ [<Px h 'i 2Ax 


The constants c* are given by 



and the smoothness measures are 


5+ =5,- [- 1 , 0 ], S+ = $[ 0 , 1 ], S^j = Sj [—2, — 1 ] , S^ = Sj[- 1 , 0 ] 


Here, 


Sj[r,s ] := Ax Y 


( A + <Pj+i\ 


»-=*■ N ' v 


^A A ( A ± <pj := ±{<fij±i - (pj). (2.27) 
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A fifth-order WENO reconstruction is obtained when d = 3. In this case, 

, + s + 2ip j+ i , 2ifj - 3 - 9<£i -2 + - 11 <Pj 

K: - 1a5 ’ W*hj~ 6Ax 

z +N _ — 2<£j-i - 3<p, + 6^j+i - ¥>j+2 , _ ~yj-2 + 6y,-_i - 3y?j - 2y J+ i 

WWaj - 6Az 

, +x + 18^+i - 9^j+2 + 2^+3 , _ 2y>j-i 4- 3y>j - 6 <pj + i + y ,+2 

C^x ) 3 ,j — >3,i - g Aa . 

The constants cj are given by 



and the smoothness measures are 

5+ = S, [-2, 0] , 5+ = S,- [-1, 1] , S 3 + = 5,- [0, 2] , 

= [ 3, —1] , [—2, 0] , 5 3 j = Sj [—1, 1] . 

The time evolution of (2.18) should be performed with an ODE solver whose order of accuracy is 
compatible with the spatial order of the scheme. In our numerical examples, we use the strong stability 
preserving (SSP) Runge-Kutta methods from [6]. 

3 Multidimensional Schemes 

In this section, we derive the two-dimensional (2-D) generalization of the compressed semi-discrete 
central- upwind scheme (2.17)— (2.18), and then extend it to three space dimensions. We also comment 
on the multidimensional interpolants that these extensions require. 

3.1 A Two-Dimensional Scheme 

We consider the 2-D Hamilton- Jacobi equation, 

<Pt + H(<p x ,ip y ) = 0, (3-1) 

and proceed as in [10]. We assume that at time t = t n the approximate point values <p* k ~ (p(xj,y k , t n ) 
are given, and construct a 2-D continuous piecewise-quadratic interpolant, <^(x,y,t n ), defined on the 
cells Sjk := {(x,y) : < 1}. On each cell S jk there will be four such interpolants (labeled 

NW, NE, SE, and SW), one for each triangle that constitutes Sj k (see Figure 3.1). Specific examples 
of <p(x , y, t n ) are discussed in §3.3. 

Similarly to the 1-D case, we use the maximal values of the one-sided local speeds of propagation 
in the x- and y-directions to estimate the widths of the local Riemann fans. These values at any grid 


point (Xj, Vk) are given by 


af k :=max{H v (ip x (x,y,t),!p y (x,y,t)) } + , 

a jk := mn{H u (& x (x,y,t),<p y (x,y,t))j_ , 

(3.2) 

bt k := max {h v ($ x (x, y,t ), tp y (x, y, f))} + , 

b~ k := min {h v ($ x (x, y, t), ip y (x, y, f)) j |, 

where C jk := x [y fc _i,y fc+ i], (•)+ 

gradient of H. 

:= max(-,0), (•)_ := min(-,0), and ( H U ,H V ) T is the 
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(j. k+j) 



Figure 3.1: Central-upwind differencing: 2-D 


The reconstruction !p(x,y, t n ) is then evolved according to the Hamilton- Jacobi equation (3.1). Due 
to the finite speed of propagation, for sufficiently small A t, the solution of (3.1) with initial data !p is 
smooth around (x£ ± ,yj ± ) where x^ ± ~ Xj ± a^At, y k± := y k ± 6^. At, see Figure 3.1. We denote by 
<Pj± f k± := Vh±> ^ n )> an d use ^e Taylor expansion to calculate the intermediate values at the next 

time level t — t n+1 


- A* • H(vJx] ±) yl ± , n, Vy(x? ± ,y n k± , t n )) + 0(Atf. 


(3-3) 


We now project the intermediate values onto the original grid points (xj,y k )' First, similarly 

to (2.5), we use new 1-D quadratic interpolants in the y-direction, ^(x™ ± , *,t n+1 ), to obtain 


^ n+1 ) = ^?±j fc- + j \+ T , - bj k - ^{^ yy )^±\ b ^ k b ^ k { At ) 2 
°jk + °jk z 


bl 


bj u 


_ b % + b tk +V*^’* + 2^ yy) ^ kb ^ k(At)2> 


(3.4) 


where ($ yy )7±\ ~ ^ys/( x "±>y]fc> * n+1 ) and Vk ; = (j/j?+ + Vk-)/ 2 - Next > we use the values 2/fc. *" +1 ) 

to construct another 1-D quadratic interpolant </>(-, t n+1 )> this time in the ^-direction, whose values 

at the original grid points are 


•p” t 1 := <t>(xj,yk,t n+1 ) 


.n+l\ 


(3.5) 


at 


= - + ^ -z-'*P(xj-,yk,t n+1 ) + + _ i/>(xj+,yk,t n+1 ) - x(<Pxx)2k la 7k a jk( At ) 2 - 


+n+l\ 


a- 


'jk 


,n *n+l\ 


'i n + 1 , 1 + /A+^2 


a jk + a jk 


a jk + a ]k 


Here, ({ p xx w y> xx (Xj,yk,t n+1 ) and := (x” + + x”_)/2. We choose (^n)^ 1 to be the weighted 
average 


( ( Pxx)j l k — , + i , - ( ( Pxx)j7- + ,+ ( (V 7 «e)j,t+i (^xx^fci ~ Vxx(%j j J/fc±! t n+1 ). (3-6) 

°jfc + °jk °jk + Ojk 


Notice that both {VyyYji^k * n (3-4) and (ifixx)]7± in (3.6) are yet to be determined. 
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We then substitute (3.4) and (3.6) into (3.5), and obtain 


n+X = “l&kVltM ± 

jk ( a jk + a jk)( b jk b jk) 


~ + t - + bj k (< pxx ) j £+ 


b jk + °jk 

btb 7, 


(At ) 2 

2 

(At) 2 


- ^ + °( At ) 2 - 


Substituting (3.3) into (3.7) yields 

+ i a 1 k bj k 

Vjk 


(4k 

+ 

a jk)(4k 

+ 

44 



a jk b jk 


( 

(4k 

+ 

a jk)(^fk 

+ 

b Jk) { 



4~k b jk 



(4k 

T 

4k)(4k 

+ 

b Jk) [ 



4k4k 



(4k 

+ 

a jk)(4k 

+ 

4k) [ 


“ j~§r [$(&,)#! + 4;.(a.)5i] ^ 

fA: ‘ °jk 

- ^ + 0(A!) 2 . 


(3.7) 


b tk b ]k 

4 k + 0 jk 

The values j. ± are computed by the Taylor expansions: 

<Pj±,k± = ± AiaJ^ ± Atb fkV$ + O(At) 2 , 

where <p*r := (p x (xj± 0, y*, t n ) and (p± := <p y ( x j> 2 /fcd: 0 , £ n ) are the corresponding right and left derivatives 
of the continuous piecewise quadratic reconstruction at (xj,yk)- 
Next, substituting (3.9) into (3.8) gives 


(3.8) 


(3.9) 


rfk 1 = <Pjk + A i4 -, jt r (y>x ~<Px) + At-j 4~r(^t ~V y )- , + ' 

3 4k + a ik ■ b jk + b jk ( a jk + a jk)( b jk + b jk) 

a jk b jk bb ( , 'P x ( X l+ > y'k+' t")’ , 'Py(~ C 'j+’ ^k+i t”)) b a jk4k^^''f ,x ^ X ]+ ) Vk—i b )> *Py ( x j+i Vk-i b )) 

+ 4k b Jk H (V*( X j-’ 2/fc+> t n ), Vy( xT j-,yk + , t")) + a+ k b+ k H (fi x (x?_^_,t n ),ipy(x]_,y%_,t n 


a % a jk 


4 k + 




(At) 2 


- k(? w )“tl + <*J,(? W )“S] ^ + 0{Mf. 


a jk + a jk 


(3.10) 
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Finally, the limit At > 0 generates a family of 2-D semi- discrete central-upwind schemes'. 

*(t) = 1 jf£l + Ml yp + a P7k H to> fy) + o t k b+ff(<p-, ip-) 

dt 3 ( a tk + a Jk)( b tk+ b Jk) 


ala, 


, J^ik_( tn + _ _\ , b jk b jk ( + _\ 

ai + a7 t ( V * ^) + 6++6-K 


U'fc T u jk 
° 7k^jk 
b jk + b jk 

b Pjk 

a jk + a jk 


'jk T 'Ufc 


bt 


^ Jk, + ^to o {A t(?ro )“«} 


(3-11) 


We still need to specify (tp xx )J jt± and {<PyyY!j7\- ^ they are proportional to (A((p x )) n+1 / Ax and to 
(A((py)) n+1 / Ay respectively, then 

Jim {At«w)£i} = 0, 


we 


and we obtain the original 2-D central-upwind scheme from [10]. However, similarly to the 1-D case, 
can choose ( <Pxx)jj±± an( ^ (pyy^^k ^ ° proportional to 1/At, so that the above limit will not vanish. 
For example, one can use the minmod limiter: 


A t(<p xx )]+l = 2 minmod ( — ~ y J±l > (&)?-, 1± 


a % + a jk 


a % + a jk 


At(£ w )"tl = 2 minmod f t n+1 ) - (<P y )"£ k . 


w vv)j±M 


b tk + b jk 


b % + b jk 


, (3.12) 


, (3.13) 


where (<p x )j± tk ± Vx(x/j ± ,yj/ ± ,t n+1 ) and (<Py)*±\ ± ■- t Py{x r j±,y'k±,t n+l ). The values of the derivative 
4> x in (3.12) are given by 


„n+l 


„n+l 


(3-14) 


A (V n .< n ^n+l'i _ < ^3+>k± Vj-k± 

<Px\,x jt y k± ,t ) , 

(aJ k + ajk )At 
and after using (3.3) and (3.9), we obtain 

M*j>Vk±,t n+1 ) = Ykl *">» Vv( x i+’ y_k ± • *")) - H (fefr?-. y£±> t n ), 5,(s?_, l/g ± , n) 

( a tk + a Jk) 


+ + + O(At). 

K+ a Jk) 


(3.15) 


Since the data is smooth along the line segments (x” ± , y% ± , t), t n <t< t n+1 , it is clear that 

Jim (<p x )]£ k± = <pt , Jim o (^x)”+i ± = ^, jun (^)SU = (3-16) 

Therefore using (3.12), (3.16), and (3.15), we obtain 


Jim {A .«(&*)£&} = minmod (V+ - ^ nt± , ^ nt± - , 


”jk t “j*. 


(3.17) 
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where 


inti a j^4. +_^fx _ -H((p x ,<p=) 

( a >* + a Jk) ( a jk + a 7fc) 

Likewise, using (3.13), we obtain 


minmod (<p+ - ^ nt± , < ti: - <p v ) , 


(3.18) 


(3.19) 


where 


int± ^jk^y __ ) 

^ ( 6 /fc + b jk ) + &7fc) 


(3.20) 

u jk _r u jkJ 

Finally, we substitute (3.17) and (3.19) into (3.11). The resulting 2-D compressed semi-discrete 
central-upwind scheme is 

_d . = g jkbpM p Z . O + a 7fc 6 j + fc g (y7 . y v ~) + (?« » yjT ) + te* » yy ) 

' JK J*' ' 

„int- 


+ °ifc a jfc 




% 






( a 7fc + a jk)( b % + &7*) 

-^ K - g T : , fat) 

ik V <& + <*,■* a Tk + <*,-* / 


u jk 


* jk ^ “jk 
„int+ 


■jk t “jfc 
,int+ 


b jk + b jk 


^ mi n m od 

-ifc V a ik + a nk + / 




Vy-Vy 

l4k +b jk 


4 k 


“ jk T “jk 
o+ _ /Jnt- 


'jk 

,int- 


JVy ~ Vy Vy “ W 

, • minmod — rr — * — rz rz 

aji. -f fljj. \ b^, + 6,- 


*jk ^ “ 'jk 

a Jk ' 

4k + a jk 


minmod 


jk 


u jk 

< u 


jk 

V~v 


4k + b jk 


(3.21) 


V + ^ 

Here, <p* n<:± and are given by (3.18) and (3.20), respectively; the one-sided local speeds, a^, and 

6& , are given by (3.2); and formulas for </?7 and are discussed in §3.3 below. 

Remark. In practice, for convex Hamiltonians 77 the one-sided local speeds are computed as 

4k ~ {#* (<?x > fy) , 0} , oJ fc = nun {bf x (<p±, <pf) , 0} , (3.22) 

b % = {-^2/ (v* . v>7) ’ °) > 4k - | ^ (vx > ^7) • °} |> 

where the maximum and minimum are taken over all the possible permutations of ±. 


3.2 A Three-Dimensional Scheme 

We consider the three-dimensional (3-D) Hamilton- Jacobi equation, 

< Pt + H (< fi x ,< Py ,< p Z ) = 0. 

We use the maximal values of the one-sided local speeds of propagation, afkl> b fkl> c %i’ * n x ~’ 
y- and z-directions, respectively. These values at any grid point {xj^y^zi) are given by the obvious 
generalizations of (3.2) and 

4ki := ma x<H w (<p x (x,y,z,t),ip y (x,y,z,t),<p z (x,y,z,t))} , 

J C jk i l ) + 
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C jkl 


mm 

Cjkl 


{H w (ip x (x, y, z, t), <p y (x t y, z, t),ip 2 (x, y, z, t))}_ 


where Cjki •'= [ x : ~i > £ J+ i] x [j/fc-i > J/fc+i] x i [ z i-i > z i+\\- Proceeding as in two dimensions, the 3-D 
compressed semi-discrete central-upwind scheme is (suppressing the indices j, k, l ) 


¥ - - p ? ^WTrW —) g 


+ ° a .- (ri -v,) + Jr^rp Wi-v,) + 


CL^~ -f- CL 


c + c 
c + + c 


- (<pT - Vz ) 


i 


(a+ + a-)(b+ + b-)(c- + c-) l'“ +< ‘ 5 [ b±c±(D i fl A,lr] 

+c + c-X) }. 


(3.23) 


where the summations are taken over all possible permutations of -f and — . For example, in the first 
sum a + b~c+ should be multiplied by H((p~ ,q>y , <pj). In (3.23), we use the notation 

: = minmod (<p+ - (<pl n \ k ±,i± > ~Vx)> 

(^>)J±h± := minmod {*t ~ " *v) ’ 

( £) ^)”±,L,; := minmod ~ (^ nt )i±,fc±,i- (^ nt )j±,jt±,i - Vz) , 


where 

/ int\ * + <Pj + Q ~^x _ H(ipZ, 

vPx )j,k±,i± • ( a + 4. a -) (a+ + a~) 

/ jnt\ = 6 + <p+ + fr<^ H(<p±,<p+,<pf) -H(ip±,ip- ,<pf) 

[(fiy W,<± : (6+ + 6-) (6+ + 6") 

/ int\ c + y+ + c~yj _ ,<pj) - ,1* ,¥>J) 

’ (c + + C“) (c + + C“) 


3.3 Multidimensional Interpolants 

The schemes developed in §3. 1-3.2 require a multidimensional non-oscillatory reconstruction. The 
simplest option is to use straightforward multidimensional extensions of the 1-D interpolants from §2.2- 
2,3, obtained via a “dimension-by-dimension” approach. 

For example, a 2-D non-oscillatory second-order central-upwind scheme is given by (3.21) with 


<Px = 


(A — i±ih T ^ , 


Ax 


Vy 


{A(f) kk±k Ay 

Ay T 2 { Vyy)j,k+ 


(<fxz)j+ i >Jb = minmod f 0 


< A ^ + §.* - (A <i,* < A <§,.-< A ^?-},. 

, ty 


(Ax) 2 


2(Ax) 2 


(Ax) 2 
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(< Pyy ) j t k + 4 = minmod 


^ (A y )" w| - (A^ w . ( A^j - (A^-j ' 


(Ax) 2 


2(Ax) 2 


(Ax) 


where 6 G [1,2], and the minmod function is given by (2.12). Similarly, the corresponding “dimension- 
by-dimension” 2-D extensions of the WENO interpolants from §2.3 can be used to reconstruct the 
derivatives in (3.18), (3.20), and (3.21). For more details see [4]. 


4 Numerical Examples 

In this section, we test the performance of the new compressed semi-discrete central-upwind schemes 
on a variety of numerical examples. We compare the methods developed in this paper, labeled BKLP, 
with the second-order scheme from [10] and the fifth-order scheme from [4], both of which are referred 
to as KNP. Our results demonstrate that the BKLP schemes achieve a better resolution of singularities 
in comparison with the corresponding KNP schemes. 

Note that in regions where the solution is sufficiently smooth, a + a~ and b + b~ are either equal to zero 
or very small (for smooth Hamiltonians and sufficiently small Ax and Ay). Hence, the BKLP and KNP 
schemes of the same order will be almost identical in these areas, and thus there will be practically no 
difference in the resolution of smooth solutions. We therefore only examine results after the formation 
of singularities, for which a + a“ and/or b + b~ may be large. 

4.1 One-Dimensional Problems 
A Convex Hamiltonian 

We first test the performance of our schemes for the Hamilton-Jacobi equation with a convex Hamilto- 
nian: 

Vt H- - {<Px + l) 2 = 0, (4.1) 

subject to the periodic initial data < p(x, 0) = — cos(7rx). The change of variables u (x, t) = <p x (x, t) + 1 
transforms the equation into the Burgers equation u t + \ (u 2 ) x = 0, which can be easily solved via the 
method of characteristics. The solution develops a singularity in the form of a discontinuous derivative 
at time t = 1/n 2 . 

The computed solutions at T = 2.5/7T 2 (after the singularity formation) are shown in Figure 5.1, where 
the second- and fifth-order BKLP and KNP schemes are compared. There is significant improvement in 
the resolution of the singularity for the BKLP schemes compared with the KNP schemes. The second- 
order BKLP scheme has a smaller error at the singularity than the fifth-order KNP scheme, while the 
fifth-order BKLP scheme has the smallest error. In Table 5.1 we show the relative L l - and L°°-errors. 

A Non- Convex Hamiltonian 

In this example, we compute the solution of the 1-D Hamilton-Jacobi equation with a non-convex 
Hamiltonian: 

(ft - cos {(p x + 1) = 0, (4.2) 

subject to the periodic initial data <p(x,0) = — cos(7rx). This initial- value problem has a smooth 
solution for t < 1.049/7T 2 , after which a singularity forms. A second singularity forms at t « 1.29/7 r 2 . 
The solutions at time T = 2/7T 2 , computed with N — 100, are shown in Figure 5.2, with a close-up 
of the singularities in Figure 5.3. The convergence results before and after the singularity formation 
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oiOiS»Ljp: conv*x •x«mpi«, N-lOO at t— 2 . 5/*x 2 



Figure 5.1: Problem (4.1). Left: The solution. Right : A close-up of the solution near the singularity. "x": 
2nd-order KNP, V’: 2nd-order BKLP, 5th-order KNP, V: 5th-order BKLP, exact solution. 



Convex example + \ (p>x + l) 2 = 0 

N 

2nd-order after singularity T = 2.5/7 r 2 

relative L l -error 

relative L°° -error 

KNP 

BKLP 

KNP 

BKLP 

100 

3.43 x io~ 4 

2.94xl0~ 4 

1.76xl0" 4 

1.29xl0 -4 

200 

4.57xl0 -5 

4.10xKr 5 

7.00 xl0~ 6 

1.96x10“° 

400 

2.15xl0“ 5 

1.85x 10 -5 

1.18xl0 -5 

8.85x10“° 

800 

2.87xl0' 6 

2.55 xl0 _e 

4.38 xl0~ 7 

1.47xl0“ 7 

N 

5th- order after singularity T — 2.5 /n 2 

relative L L -error 

relative L°° -error 

KNP 

BKLP 

KNP 

BKLP 

100 

1.50xl0~ 4 

1.13 xlO -4 

1.46 xlO" 4 

l.lOxlO -4 

200 

2.33x10-° 

9.95xl0 -7 

1.56xlO _B 

3.42xl0“ ? 

400 

9.39x 10 - ° 

7.08 x 10~ 6 

9.33xl0 _e 

7.02x10“° 

800 

1.13xl0~ 7 

3.94xl0 -8 

7.54xl0 -s 

1.41x10“° 


Table 5.1: Problem (4.1). Relative L 1 - and T°°-errors for the KNP and BKLP schemes. 
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ikxwwivw .xampta, N-100 » 



Figure 5.2: Problem (4.2) — the KNP and BKLP numerical solutions. 

closeup: non>conv«x •xampla, N^IOO to 1<»2.0 /k 2 clo*«up: non-conv#x N«lOO to t*a.O/ji 2 




Figure 5.3: Problem (4.2). Left: The singularity near x = 0.25. Right: The singularity near x = 1.11. "x": 
2nd-order KNP, “o” : 2nd-order BKLP, : 5th-order KNP, V: 5th-order BKLP, exact solution. 


axe given in Table 5.2. In this example, the local speeds of propagation were estimated by (2.2). The 
results are similar to the convex case, though the improvement here is somewhat less dramatic. 

Next, we examine the convergence of the numerical solutions of (4.1) and (4.2), computed by the 
fifth-order BKLP and KNP schemes. These results, together with the fifth-order methods from [7] and 
[3], are shown in Figure 5.4. The reader may note that the convergence rates in these examples are 
erratic. However, this investigation of the relative L ^errors for many different grid spacings shows that 
the behavior is due to super-convergence at some grid spacings. Notice that for all grid spacings the 
L x -e rror of the BKLP method is less than the others. 


4.2 Two-Dimensional Problems 

In this section, we test the 2-D BKLP schemes on Hamilton-Jacobi equations with convex and non- 
convex Hamiltonians. We start with the convex problem (compare with (4.1)): 


<Pt + g (*Px + Vy + 1) = 0> 


(4.3) 





relative L -error 
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Non-convex example ip t - cos (tp x + 1) = 0 

N 

2nd- order after singularity T — 2.0 /n 2 

relative L l -error 

relative L 00 -error 

KNP 

BKLP 

KNP 

BKLP 

100 

5.59xl0 -4 

4.97xl0 -4 

1.75 x 10~ 4 

1.22xl0" 4 

200 

9.52X10" 5 

9.52X10" 5 

4.47x10"° 

4.11x10"° 

400 

2.40xl0 -5 

2.40xl0 -5 

2.06x10"° 

1.57x10"° 

800 

6.02xl0 -6 

6.02x10"° 

6.30xl0" 7 

3.09xl0" 7 

N 

5th-order after singularity T = 2.0/7 r 2 

relative L l -error 

relative L°°- error 

| KNP 

BKLP 

KNP 

BKLP 

100 1 

1.48xl0" 4 

9.91xl0" 5 

1.25X10" 4 

8.17xl0" 5 

200 

8.49x 10 -8 

5.82xl0" s 

1.35X10" 7 

8.88x10"° 

400 

7.89x 10 -9 

6.60xl0" 9 

8.21 xlO" 7 

5.48xl0" 7 

800 1 

6.63x 10 -1 ° 

5.13xl0~ lu 

2.25xl0" 7 

7.77x10"° 


Table 5.2: Problem (4.2). Relative L 1 - and T°°-errors for the KNP and BKLP schemes. 


convex H, T— 1 ,5/rc 2 




number of points 


Figure 5.4: Convergence of the 1-D examples. Left: Problem (4.1), T = 2.5/tt 2 . Right: Problem (4.1), 
T = 2/? r 2 . 5th-order BKLP, “A”: 5th-order KNP, “o": the 5th-order method from [7], The solid lines 

show example rates of convergence. 
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which can be reduced to a 1-D problem via the coordinate transformation 

( = ( I/ 2 1/2 \f x \ 

W V I/ 2 1/2 J\yJ 

The relative Z 1 - and Z°°-errors for the periodic initial data p (x, y, 0) = — cos ( 7 r(x + y)f 2) = — cos (7 t£) 
after the singularity formation at T = 2.5/7 r 2 are shown in Table 5.3. 



2-D Convex example 

N 

2nd-order after singularity T = 2.5/7T 2 

relative L l -error 

relative Z°° -error 

KNP 

BKLP 

KNP 

BKLP 

50 

7.31 xl0~ 4 

7.26x10^ 

2.13xl0“ 6 

2.23xl0' fa ' 

100 

3.27xl0~ 4 

2.99xl0 -4 " 

1.58xl0 _ti 

1.30xl0~ 6 I 

200 

4.37xl0 -5 

4.12xl0~ 5 

2.34xl0“ 8 

9.87xl0~ y 

N 

5th- order after singularity T = 2.5/7r 2 

1 relative L l -error 

relative Z°° -error || 

KNP 

BKLP 

KNP 

BKLP 

50 

6.01 xl0~ 5 

4.39xl0 -4 

3.79x10^ 

1.60xl0~ 7 

100 

1.40xl0- 4 

1.19xl0 -4 

1. 33x10“° 

l.llx 10~ 6 

200 

1.98xl0 _ ° 

1.23x10-° 

5.97xl0~ y 

2.58xl0~ y j 


Table 5.3: Problem (4.3). Relative Z 1 - and Z^-errors for the 2-D KNP and BKLP schemes. 


In Table 5.4, we present similar results for the non-convex problem (compare with (4.2)): 

(ft - COS (p x + (fy + 1) = 0, 

with the periodic initial data p (x, y, 0) = — cos ( 7 r(x + y)/2). 


(4.4) 



2-D non-convex example 

N 

2nd- order after singularity T = 2.0/7 r 2 

relative Z 1 -error 

relative Z°° -error 

KNP 

BKLP 

KNP 

BKLP 

50 

1.84xl0~ 3 

1.75xl0“ 3 

5.11xl0 _6 

3.66 xl0~ 6 

100 

5.86xl0“ 4 

5.48xl0“ 4 

1.50x10-° 

1.21x10-° 

200 

1.14xl0“ 4 

1.13xl0~ 4 

2.15x10-° 

2.04x10-° 

N 

5th-order after singularity T = 2.0/7 r 2 

relative L l -error 

relative Z°°- error 

KNP 

BKLP 

KNP 

BKLP 

50 

1.79xl0“ 4 

1.44xl0“ 4 

6.81xl0- r 

6.80xl0- 7 

100 

1.14xl0 -4 

8.97xl0“ 5 

1.02x10-° 

7.82 xl0~ 7 

200 

4.42 xlO- 7 

4.32xl0 -7 

5.65xl0 _lu 

4.18x10-™ 


Table 5.4: Problem (4.4). Relative Z 1 - and Z°°-errors for the 2-D KNP and BKLP schemes. 
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A Appendix: A Proof of Theorem 2.1 


Proof. First, we fix u and show that H BKLP (u,u ) given by (2.22) is a non-increasing function of u 
(the proof that H BKLP is non-decreasing function of its second argument is similar). We denote by 


Q(u) := { 


H(u) - H(u ~ ) 


u — u 




, U^U , 


u = u , 


(A.1) 


and 


1 


A{u) := - (a + (u) - a (u)) , 


where a + {u) = ),0} and a (u) = | min^^u), H'iu ),0}|, and by U\, U 2 , Vi, V 2 , 

the sets 

U x := Ui{vT) = {u : Q(u ) - A(u ) < 0}, U 2 := U 2 (u^) = {u : <9(u) - A(u) > 0}, (A.2) 

V x := T^(u-) - {u : Q(u) - A(u) < 0}, V 2 := V 2 (u~) = {u : Q(u) - A(u) > 0}. (A.3) 

Both U x and V 2 are open sets (Q and A are continuous) and as such can be represented as a union of 
at most countably many disjoint open intervals Ij and J}, respectively, i.e. 


U 1 = U JL x Ij and V 2 = J 3 . 

In the new notation it is easy to verify that the Hamiltonian H BKLP can be written as 

H bklp {u ) u~) y u e 

Hf KLP (u,u~), ueU 2 , 


H bklp (u,u ) := 


or as 


H bklp (u,u~ ) := 


H bklp (u, U ), u € Vi, 

H? klp (u,u~), uev 2 , 


(A.4) 


(A.5) 


(A.6) 


where 


zsBKLPr 0 {u)H(u + )+a + {u)H{u ) a + (u)a (tt) 
n i l — ztt77\ , rr7T\ i+rrrr /..vw — u ) 


a + (u) + a~(u) 
a + (u)a~(u) 


(a+(u) + a~(u)) 2 


a + (u) + a~(u) 
(u - u~) [Q(u) + a - (u)] , 


(A.7) 


and 


a + (u) + a~(u) a + (u) + a~(u) K 


+ 


a + (u)a (u) 
(a + (u) + a~(u))‘ 


:(u — u ) [a + (u) — Q(u)] . 


(A.8) 


Notice that for those u for which Q(u) = A(u), we have that H BKLP (u,u ) = H BKLP (u,u ), and 
therefore H BKLP {-, u~) is a well defined function consisting of the pieces H BKLP (-, u~), i — 1,2. We will 
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use formulas (A.5) and (A. 6) and continuity arguments to show that H BKLP (-, u~) is a non-increasing 
function on the whole real line. 

Consider the point u * such that H'(u*) — 0 (see assumption (Al)). Then H BKLP {u,u~') ) i = 
1,2, are continuously differentiable on the intervals (— oo, min(u _ ,u*)), (min(u _ ,u*),majc(u _ ,u*)), and 
(max(u“, u*), oo) (in case u~ — u*, on the first and third interval only) and continuous on JR. 

Case 1. Let u G (— oo,min (u',^)). Then ^ (a + (u)) = 0 since a + (u) = max{H'(u) y H'(u~),0} and 
H f is a non-decreasing function of u. In this case we have 




2(a {u)) l a (u)a + (u)(u — u ) 
(a+(u) + a"(u)) 3 
(a-(u)f(a+(u)-H>(u)) 
(a + (u) + a~(u)) 2 


[a + (u) - <5(u)j 


Note that there exists £ G (u ) u ) such that 


Q(u) := 4 ( 4 = H ’^) < H '(vT) < a+(u), 

u — u~ 

a~ is a smooth non-increasing function on (— oo,min(u”,u*)), a + (u) > 0, a~~(u) > 0, H f (u) < a + (u), 
and therefore the derivative (H 2 KLP (u, u~)) < 0. Hence H BKLP {u^u~) is non-increasing on 

(— oo,min(u“,^*))- 

Similarly, for u G (— oo, min(u“, u*)) we have 


d , r bklp 
au K 


(u,u )) 


2(q-(tt)) / (q+(tt)) 2 (u-u-) 
(a + (u) + a~(u)) 3 
a~(u)a + (u)(a+(u) — 
(a+(u) + a“(w )) 2 


[Q(u) - A{u)\ 

g-(u)H'(u) 
a+(u) 4* a~(u) ’ 


(A.9) 


Now we fix j and consider the corresponding open interval Ij n (— oo, min(u~, u*)) (see (A.4) and 
representation (A.5)). As above, a~ is a smooth non-increasing function on (— oo, min(u”, u*)), a + (u) > 
0, a~{u) > 0. On each Ij , Q(u) < A(u), and therefore the first term on the right-hand side (RHS) 
of (A.9) < 0. The second term is non-positive since a + (u) > H f (u). The last term is < 0 because 
H'(u) < 0 for u e (—'oo, mi n(u",u*)). This proves that H PKLP (u,u~) is a non-increasing function of 
u on Ij n (— oo,min(u”,n*)), for every j. 

Case 2. Let u G (rnin(u“,u*), max(u“,u*)). In this case the derivatives are ^(a - ^)) = 0 and 
^ (a~{u)) — 0. Therefore 


I- ( Hi KLP (u,u “)) = -HW < 0 

du v ' (a + (u) + a~ (u)) 2 


since a + (u ) > H'(u), 


which shows that H PKLP {-,u ) is non-increasing on (min(u ,tt*),max(u Likewise 

A (tjBKLP ( . -s\ fl~(n)a + (u)(a + (n) a~ («)#'(«) 

du ' 1 ,U (a + (u) + a~ (u)) 2 a + (u) + a~(u) 

The first term on the RHS is < 0 since q + (u) > H'(u). As for the second term, we have two possibilities. 
If u~ < u < u* then H'(u) < 0, which will make the whole term <0. If u* < u < u~, then a~ (u) = 0 
and the second term is 0. Therefore, the derivative of H PKLP is < 0, and thus H PKLP (-, u~) is non- 
increasing on (min(u~,u*),max.(u~,u*)), and in particular on Ij fl (min(u - , u*), max(n _ , u*)) for every 
J f- 
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Case 3a. Let u* < u~ < u. Then a~(u) = 0 and therefore H BKLP = H BKLP (u,u~) = 
H BKLP (u,u~) = H(u~). In particular, H BKLP is a non- increasing function on (u“,oo), and H BKLP 
is a non- increasing function on Ij H (ti“,oo), for every j. 


Combining the results from Cases 1, 2 and 3a. we obtain that H BKLP is non-increasing on the whole 
real line (since it is continuous on JR and non-increasing on each of the intervals (— oo, mir^u - , u*)), 
(min(u~,u*), max(u“, u*)), and (max(u”,u*),oo)), and H BKLP is a non-increasing function on every 
open interval Ij from U\ (same reasoning). Since H BKLP (u,u~) = H BKLP (u,u~) for u e dlj it will 
follow from (A. 5) that H BKLP (' , u ) is non-increasing on the whole real line. 

Case 3b. Let u~ <u* < u . In this case we will utilize representation (A. 6) for H BKLP and, using 
the results from Cases 1 and 2, we will show that H BKLP is a non-increasing function on the interval 
[a, b ] for any a and 6, and therefore on the whole real line. 

Notice that in this case Vi ~ S~(u,u~), and then, by assumption (AS), V\ fi [a, b] is either 0 or a 
finite number of points and/or a finite union of closed intervals 7^. Note also that we have 


( u *, oo) f) [a, b] = [Jj ft ( u * , oo) D [a, 6]] U U ™ = {Tk, for some m. (A. 10) 

For u e (u*,oo) we have that J- (a~(u)) = 0, and hence 

i- (Hr LF («,«-)) = _ AM] _ (''{ffw-y . 

du v 7 (a+(u) + u~(u)) 3 (a + (w) + a~(u)) 2 

As in Case 1, we fix j and consider this time the corresponding interval Jj Pi (u*, oo). Since a + is a 
smooth non-decreasing function on (u*, oo) and Q(u) > A(u) on Jj , the first term on the RHS < 0. Also 
a+(u) > H f (u ) and hence the second term is also < 0. This gives that H BKLP (u , u~) is a non-increasing 
function of u on Jj fi (u*, oo) for every j. 

When vT < u* < u, a + (u) = H f (u ), a~(u) = - H f (u ~ ), and hence 




a (u)H'(u) 
(a+(u) + a _ (u )) 3 


G(u,u ), 


(A.11) 


where G(u,u~) is given by (2.20). Since H f (u) > 0 for u > u *, conditions (2.20)-(2.21) ensure that 
the RHS in (A.ll) is < 0 for u G 7U This shows that H BKLP (u,u~) is a non-increasing function on 
each of the intervals 7& constituting Vi 0 [a, b] (if Vi fl [a, 6] consists of finite number of points, then 
ffBKLP = h bklp at these points). 

Since H BKLP (u , u~) = H BKLP (u , u~) on dJj , all the above arguments and (A.6) prove that H BKLP 
is non-increasing on (u*, oo). This, together with the conclusion in Cases 1 and 2 and the continuity of 
HP KLP (u,u~), i~ 1,2, give that H BKLP is non-increasing on [a, b]. 


Similarly, one proves that H BKLP (u+ , u) (when u + is fixed) is a non-decreasing function of the second 
argument u. Here, in the case corresponding to (A.ll) in Case 3b above, we have 




«)) = 


a + (u)H'(u) 
(a+(u) + a~(u )) 3 


G(u,u + ), 


u < u* < u + . 


Since H'(u ) < 0 for u < u*, conditions (2.20)-(2.21) guarantee that the derivative is non-negative, and 
hence H£ klp (u + ,u) is a non-decreasing function of u on the finite union of closed intervals S + (u, u + ) fl 
[a,b]. 
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